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1.1 


DISCLAIMER 

The  reader  is  cautioned  that  computer  programs  developed  in  this  research 
may  not  have  been  exercised  for  all  cases  of  interest.  While  every  effort  has 
been  made,  within  the  time  available,  to  ensure  that  the  programs  are  free  of 
computational  and  logical  errors,  they  cannot  be  considered  validated.  Any 
application  of  these  programs  without  additional  verification  is  at  the  risk  of  the 
user. 
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I.   INTRODUCTION 

In  modern  missile  warfare,  new  technologies  are  producing  faster,  more 
accurate  anti-ship  cruise  missiles  (ASCM).  They  are  able  to  fly  at  very  low 
altitudes  and  perform  a  high  speed  'pop-up'  maneuver  just  before  impact,  in 
order  to  defeat  the  ship's  missile  defenses,  as  shown  in  Figure  1.  In  order  to 
counter  the  advanced  ASCM's,  new  control  schemes  must  be  developed  to 


Figure  1.  Anti-ship  cruise  missile  terminal  maneuver 

enable  the  ship's  surface-to-air  missiles  to  counter  this  threat. 

The  aim  of  this  thesis  is  to  illustrate  the  concept  of  Bang-Bang  control  and 
how  it  applies  to  control  of  vertically-launched  surface-to-air  missiles.  Bang- 
Bang  control  can  be  used  in  the  missile  boost  phase  in  order  to  turn  the  missile 
downrange  in  minimum  time,  and  also  in  the  terminal  phase,  where  the  missile 
must  react  very  quickly  in  order  to  intercept  a  maneuvering  target. 

This  work  applies  the  minimum-time  (Bang-Bang)  control  concept  to 
vertically-launched  surface-to-air  missiles.     In  Chapter  II,  the  concept  of 


control  is  defined  and  solved  for  a  second  order  system.  In  Chapter  III,  the 
equations  of  motion  for  a  vertically- launched  missile  in  the  boost  and  terminal 
phases  are  developed.  In  Chapter  IV,  the  boost  simulations  comparing  missile 
trajectory  versus  maximum  thrust  vector  control  (TVC)  angle  are  presented. 
In  Chapter  V,  the  terminal  phase  missile  simulations  are  presented  for  a  non- 
maneuvering  target  and  also  for  a  maneuvering  target. 


II.  CONCEPT  OF  BANG-BANG  CONTROL 

In  this  section  the  form  of  the  optimal  control  for  a  particular  class  of 
systems  will  be  determined  by  using  Pontryagin's  Minimum  Principle. 

A.  PONTRYAGIN'S  MINIMUM  PRINCIPLE 

It  shall  be  assumed  that  the  state  equations  of  the  system  are  of  the  form 

x  =  Ax  +  Bu,  (2-1) 

where  A  is  an  n  by  n  array,  B  is  an  n  by  m  array,  and  A  and  B  may  be 
explicitly  dependent  on  the  states  and  time.  For  this  research  it  will  be  assumed 
that  there  is  a  single  input;  therefore,  m  is  1.  It  is  specified  that  the  admissible 
controls  must  satisfy  the  inequality  constraints 

N.  <  u(t)  <  N+  (2-2) 

where  N.  and  N+  are  the  known  lower  and  upper  bounds  for  the  control  input. 
It  is  desired  to  drive  x  from  an  initial  state  x(to)  to  to  a  desired  final  state  x(tf), 
where  to  is  the  problem  start  time  and  tf  is  the  problem  end.  The  optimal  control 
is  that  control  which  drives  the  state  from  its  initial  state  to  desired  final  state 
using  the  least  "cost".  The  cost  can  be  any  desired  measure  of  system 
performance.  In  the  case  of  a  minimum-time  controller,  the  cost  function  can 
be  represented  as 

j  =  Jdt  =  tf-t0. 

l°  (2-3) 


Pontryagin's  Minimum  Principle  states  that  the  optimal  control,  u*,  which 
minimizes  the  cost  function,  must  minimize  the  Hamiltonian,  which  is  defined  as 

H(x(t),  u(t),  p(t),  t)=  l+pT(t)(Ax(t)+Bu(t»  (2-4) 

where  p(t)  is  a  Lagrange  multiplier  vector,  and  is  arbitrary.  p(t)  can  be  defined 
as  the  'costates'  of  the  system,  and  p(t)  can  be  written  as 


P(t)  = 


Pi(t) 

i_p2(t)_r 


(2-5) 

The  necessary  condition  for  the  optimal  control  u*  to  minimize  the  cost 
function  J  is 

H(x*(t),u*(t),p*(t),t)  <  H(x*(t),u(t),p*(t),t)  (2-6) 

for  all  times  between  to  and  tf  and  for  all  admissible  controls.    The  asterisks  in 
equation  (2-6)  represent  the  optimum  values,  and  thus, 

l+p*T(t)Ax*(t)+p*T(t)Bu*(t)>l+p*T(t)Ax(*t)+p*T(t)Bu(t)  (2-7) 

Therefore 

p*T(t)Bu*(t)  <  p*T(t)Bu(t)  (2-8) 

for  all  admissible  u(t)  and  for  all  times  from  tinitial  to  t<-mal.  If  u(t)  is  constrained  to 
±N,  the  optimal  control  U*(t)  is 

-•(»)-{-"        ;°rr°U-Nsign(p;)  (2-9) 

l+N        forp2<0J  v     ' 


This  indicates  that  the  time-optimal  control  is  'Bang-Bang';  that  is,  the 
optimum  control  switches  between  its  maximum  positive  and  negative  values. 
There  are  also  three  ideas  which  deserve  mention: 

(1)  (Existence) 

If  all  the  eigenvalues  of  A  have  nonpositive  real  parts,  then  an  optimal 
control  exists  that  transfers  any  initial  state  xo  to  the  origin.  If  there 
are  positive  eigenvalues  (i.e.  unstable  roots),  there  may  be  some 
region  of  the  state  space  where  the  system  is  uncontrollable. 

(2)  (Uniqueness) 

If  an  optimal  control  exists,  then  it  is  unique. 

(3)  (Number  of  Switchings) 

If  the  eigenvalues  of  A  are  all  real,  and  a  (unique)  time-optimal 
control  exists,  then  the  control  can  switch  at  most  (n-1)  times. 
Thus,  an  nth-order  system  having  all  real,  non-positive  eigenvalues  has  a  unique 
time-optimal  control  that  switches  at  most  n-1  times  (not  counting  switching  off 
at  time  tfinai).  For  complex  conjugate  eigenvalues,  more  than  n-1  switchings 
may  be  required.[l] 

B.     APPLICATION  TO  A  SECOND  ORDER  SYSTEM 

In  order  to  illustrate  the  minimum  principle,  consider  the  second  order 
differential  equation 

x  =  u,  (2-10) 

or  in  state  space, 


LA2. 


"0     1 

0    0 


Txi]  |~o' 

x2      [l 


(2-11) 


From  the  minimum  principle,  the  time-optimal  control  for  this  system  is  ±N. 
Thus  the  segment  of  optimal  trajectories  can  be  found  by  integrating  the  state 
equations  (2-11),  with  u=±N,  from  time  to  to  t.  This  yields  the  results 


x2(t)  =  ±Nt  +  C1 


(2-12) 


X!(t)  =  ±N— +C1t+C2 


(2-13) 


where 


C1=+Nt0  +  x2(t0) 


(2-14) 


and 


C2=TN^-C1t0  +  x1(t0). 


(2-15) 


Ci  and  C2  are  functions  of  the  initial  conditions  and  can  be  treated  as  constants. 
Time  can  be  eliminated  from  equations  (2-12)  and  (2-13)  by  squaring  the  first 
equation,  dividing  the  result  by  2N,  and  comparing  the  result  with  equation 
(2-13)  to  obtain 


x1(t)  =  ^rx22(t)  +  C3  forU  =  +N 


(2-16) 


xi(t)  =  -^-x22(t)  +  C4 


forU  =  -N 


(2-17) 


where 


Co   —  Ct 


2N 


(2-18) 


C4=C2+-^- 

4       2     2N 


(2-19) 


C3  and  C4  are  functions  of  the  initial  conditions  and  may  also  be  treated  as 
constants.  Equations  (2-18)  and  (2-19)  each  define  a  family  of  parabolas,  as 
shown  in  Figures  2  and  3. 

Analysis  of  Figures  2  and  3  reveal  the  controls  corresponding  to  the 
following  situations: 

1.  u*=+l  implies  that  the  initial  state  x0  lies  on  segment  AO  at  time  to. 

2.  u*=-l  implies  that  the  initial  state  x0  lies  on  segment  BO  at  time  to. 

This  defines  the  optimal  control  to  reach  the  origin.  In  order  to  reach  these 
zero  trajectories,  the  opposite  value  of  control  is  required.   It  can  be  seen  from 
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Figure  2.     Trajectories  for  u=-N 
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Figure  3.     Trajectories  for  u=+N 

Figure  4  that  the  two  segments  AO  and  BO  form  a  single  curve.  This  curve  is 
known  as  the  switching  curve.  All  states  that  fall  above  the  switching  curves 
will  result  in  the  control  input  U=-N,  and  those  states  which  fall  below  the 
switching  curves  will  result  in  the  control  input  U=+N.  This  curve  can  be 
represented  as 


1 


X!(t)  =  -— -X2(t)|x2(t)| 


2N 


(2-20) 


Thus,  the  control  law  can  be  obtained  by  moving  all  the  terms  of  equation 
(2-20)  to  one  side: 

1 


U(t)  =  -Nsign  x1(t)  +  — -x2(t)|x2(t)| 


2N 


(2-21) 
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Figure  4.  Optimal  switching  curves 


Examples  of  the  use  of  this  switching  law  is  shown  in  Figure  5.  From  any 
arbitrary  initial  state,  the  trajectory  will  follow  the  parabola  that  passes  through 
the  start  point  until  the  trajectory  reaches  the  optimal  switching  curve.  Then  the 
control  switches  signs,  and  follows  the  optimum  parabolic  path  to  the  origin. 

Equation  (2-21)  is  the  mathematical  representation  of  the  Bang-Bang 
control  law.  This  control  law  can  be  used  in  numerous  applications  where  it  is 
desired  to  transfer  a  state  from  initial  conditions  to  desired  final  conditions  in 
minimum  time.  Two  examples  that  are  of  interest  in  this  research  are: 

1)  drive  a  state,  such  as  a  position,  to  zero  in  minimum  time.  The  control 
law  for  this  portion  of  the  missile  flight  is  of  the  form 


U(t)  =  -Nsign(^x1(t)  +  ^-x2(t)|x2(t)|\ 


(2-22) 


Ill                                         III 
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Nv          J               xi 

Figure  5.  Optimal  control  switching  schemes 

where  xi(t)  is  the  position  at  time  t,  x2(t)  is  the  velocity  at  time  t,  and  N  is  the 
maximum  admissible  force  that  can  be  used  to  accelerate  or  decelerate  the 
state;  and 

2)  drive  the  time  rate  of  change  of  a  state  to  zero.  In  this  case  it  is 
desired  to  keep  a  state  constant,  as  opposed  to  driving  it  to  zero.  The  form  of 
the  controller  is  also  given  by  equation  (2-22),  but  X!  in  this  case  will  represent 
velocity  and  x2  will  represent  acceleration,  with  N  being  the  maximum 
admissible  control  with  which  to  change  the  states. 

In  the  next  chapter,  the  application  of  the  second  order  Bang-Bang 
controller  to  boost  phase  and  terminal  phase  missile  control  will  be  examined. 
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III.     DEVELOPMENT  OF  EQUATIONS  OF  MOTION  FOR  A 
VERTICALLY-LAUNCHED    MISSILE 

The  equations  of  motion  for  a  missile  can  be  derived  from  Newton's  second 
law  of  motion,  which  states  that  the  summation  of  all  external  forces  acting  on  a 
body  must  be  equal  to  the  time  rate  of  change  of  its  momentum,  and  the 
summation  of  the  external  moments  acting  on  a  body  must  be  equal  to  the  time 
rate  of  change  of  its  moment  of  momentum  (angular  momentum).  The  time 
rates  of  change  can  be  expressed  by  two  vector  equations: 

2*-|Mt),  (3-D 


and 


YMi=i  (3-2) 

^     '      dt 


where  i  indicates  ith  cartesian  coordinate  of  the  vector  with  respect  to  inertial 
space.  By  definition,  H  is  the  angular  momentum,  or  moment  of  momentum,  of 
a  revolving  body.  [2] 

There  are  several  assumptions  that    are  made  in  order  to  simplify  the 
problem: 

1 .  The  mass  of  the  missile  remains  constant. 

2.  The  missile  airframe  is  a  rigid  body. 

3.  The  earth  is  an  inertial  reference,  and  unless  otherwise  stated,  the 
atmosphere  is  fixed  with  respect  to  earth. 
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In  order  to  completely  describe  the  motion  of  a  missile  in  three  dimensions, 
a  total  of  six  nonlinear  differential  equations  must  be  solved.  However,  certain 
assumptions  will  be  made  in  order  to  reduce  the  number  of  equations  needed  to 
adequately  describe  the  motion  of  the  missile: 

1.  The  X  axis  will  be  assumed  to  be  downrange  from  the  launch 
platform  to  the  target,  the  Y  axis  will  be  crossrange,  and  the  Z  axis  will  be 
altitude. 

2.  The  type  of  missile  studied  here  is  known  as  a  skid-to-turn  missile, 
as  this  missile  uses  direct  side  force  to  turn.  Thus,  there  will  be  no  coupling 
terms  between  the  pitch  and  yaw  axes,  and  one  set  of  equations  will  describe 
the  motion  in  each  of  the  Y  and  Z  planes. 

3.  It  will  be  assumed  that  the  missile  is  roll  stabilized,  i.e.,  that  there 
will  be  no  turning  moments  about  the  X  axis. 

4.  Since  there  are  no  cross-coupling  terms  between  the  pitch  and  yaw 
axes,  the  simulation  can  be  broken  into  two  problems,  one  where  the  pitch 
angle  is  held  constant,  and  the  other  where  the  yaw  angle  is  held  constant.  In 
this  discussion  the  yaw  angle  will  be  assumed  to  be  constant  at  zero.  Using 
these  assumptions,  a  3  degree-of-freedom  model  can  be  developed.  The  force 
diagram  of  a  3  DOF  missile  in  flight  is  shown  in  Figure  6.  The  governing 
equations  of  motion  of  this  missile  are: 

XFx=mx=  -Fliftsinym  +  Tcos(0m+5)  -  Fdragcosym  (3-3) 

£Fz=mz  =  -mg  +  F^cos^  +  Tsin(Gm  +  5)  -  Fdragsinym        (3-4) 
£Mcg  =  icge  =  -FliftLpgcosym-FdragLpgsinYm-TLtgsin5         (35) 
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Lift                                                   /*\ 

\                 xx    y\a        v      \ 

Center  of           \                /  /  s^^00**                \ 
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Figure  6.    Forces  acting  on  a  missile  in  flight 

m  missile  mass  (Kg) 

g  gravity  (9.81  m/s^) 

Fiift  aerodynamic  lift  force  (N) 

Fdrag  aerodynamic  drag  force  (N) 

Ym  angle  of  velocity  vector  from  reference 

0m  angle  of  missile  head  from  reference 

5  TVC  deflection  angle 

leg  moment  of  inertia  of  the  missile  about  the  center  of  gravity 
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Lpg      distance  from  aerodynamic  center  of  pressure  to  center  of 

gravity 
Ltg       distance  from  TVC  control  to  center  of  gravity 
The  aerodynamic  force  acting  on  a  lifting  body,  in  this  case  the  wings  and 
control  surfaces  of  the  missile,  may  be  resolved  into  two  components,  which 
are 

Fiift=CL|v*S  (3-6) 

and 

Fdrag=CD-V^S  (3?) 

where  (p/2)Vm2  is  the  dynamic  pressure  in  N/m2. 

The  lift  and  drag  coefficients  are  functions  of  angle  of  attack  (a),  and  Mach 
number.  Actual  values  for  Cl  and  Co  must  be  measured,  as  they  are  highly 
dependent  on  the  airfoil  geometry.  However,  Cl  and  Cd  may  be  modeled 
satisfactorily  using  the  following  approximations: 

CL=0.1a  (3-8) 


C^^h(l  +  0-2a2)-  0-9) 

r     max 


These  coefficients  are  in  reasonable  agreement  with  reference  [2]  at  low 
angles  of  attack  (a<10°).  The  only  time  the  missile  can  expect  to  exceed  10° 
angle  of  attack  is  during  the  tip-over  phase,  when  velocity,  and  therefore  lift  and 
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drag  forces,  are  small  compared  to  the  thrust  force.  This  is  considered  an 
appropriate  approximation  for  the  scope  of  the  model  under  discussion. 

A.  BOOST  PHASE  CONTROLLER 

The  primary  objective  of  the  missile  control  system  in  the  boost  phase  of 
flight  is  to  tip  the  missile  over  and  point  the  missile  in  the  general  direction  of  the 
target  in  minimum  time.  Therefore,  6,  the  angle  that  the  missile  makes  with  the 
horizontal  reference,  is  a  logical  choice  for  a  control  variable.  If  the  force  that 
creates  the  moment  that  tips  the  missile  over  can  be  controlled,  a  general 
control  algorithm  can  be  formulated.  This  is  shown  graphically  in  Figure  6. 

At  launch,  the  missile  is  assumed  to  be  in  a  vertical  position,  with  the 
velocity  vector  pointing  straight  up,  i.e.  6  =  90°  and  y  =  90°.  In  the  initial  boost 
phase,  it  is  assumed  that  the  first  two  terms  of  equation  (3-5)  are  negligible  as 
compared  to  the  third.  Therefore,  equation  (3-5)  may  be  simplified  to 

-TL„„sin8 


6  = 


-eg 


I 


eg 


(3-10) 


In  state  space,  this  can  be  represented  as 


0     1 
0    0 


V 


TLcgsin8 


A 


I 


eg 


(3-11) 


The  only  term  in  equation  (3-11)  that  can  be  used  to  control  9  is  5,  the  TVC 
deflection  angle.  A  simple  and  effective  method  for  controlling  missile  tip  over 
is  to  employ  the  Bang-Bang  controller,  developed  in  Chapter  II,  to  drive  0  to  0 
(horizontal  flight)  in  minimum  time,  as  shown  in  Figure  7.  Thus,  the  control  law 
is 
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V 


\ 


Center  of 
Gravity 


TVC  Force  Perpendicular 
to  Missile  Body 


Figure  7.  Boost  Phase  Diagram 


u(t)  =  -Nsignfe(t)+^~e(t)|e(t) 


(3-12) 


where  N  is  the  maximum  torque  available  to  rotate  the  missile  and  is 
determined  from  the  relation 


N  =  -^-sin(5max) 


'Cg 


(3-13) 


and  5max  is  the  maximum  thrust  vector  control  (TVC)  angle  available. 

This  indicates  that  for   minimum  tip-over  time,  Smax=90°,    However,    this 
may  not  be  desirable  since  all  of  the  thrust  will  go  to  rotate  the  missile,  while 
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none  of  the  thrust  will  accelerate  the  missile  downrange.  This  amounts  to  the 
missile  'pin-wheeling'  in  the  air.  Therefore,  TVC  angles  of  less  then  90°  should 
be  used.  Several  TVC  angles  will  be  examined  to  show  the  tradeoffs  involved. 
When  the  missile  achieves  horizontal  orientation  (i.e.  0=0°),  the  control  law 
will  continue  to  switch  from  +N  to  -N  in  an  attempt  to  keep  the  missile  at 
exactly  0°.  This  phenomenon  is  known  as  'chattering'.  At  this  point  it  is  no 
longer  desirable  to  use  maximum  control  effort  to  correct  for  small  changes  in 
9.  A  satisfactory  method  for  removing  the  chatter  from  a  Bang-Bang 
controller  is  to  employ  linear  control  in  a  small  region,  of  width  e,  about  the 
desired  value,  as  shown  in  Figure  8. 
This  adapted  Bang-Bang  controller  can  be  expressed  mathematically  as 


u(t)  = 


Nsignfe(t)+^e(t)|e(t)|l      e(t)+^_e(t)|e(t) 


>£ 


—  e(t)+— -0(t)e(t) 

e  V  2N  J 


-e<e(t)+^-e(t)|8(t) 


<£ 


(3-14) 


1 

+N- 
-e 

+u(t) 

X, 

»»  ii(n\ 

\    v 

1 

+e 

— N 

—  ■  -^  1  J[Q) 

\ 

> 

Figure  8.  Bang-Bang  Control  Law  with  Linear  Zone 
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B.  TERMINAL  PHASE 

1.   Missile-Target  Geometry 

In  this  section,  it  will  be  assumed  that  the  vertically  launched  surface- 
to-air  missile  has  completed  tip-over,  and  has  settled  out  to  a  constant  speed 
and  altitude.  The  start  point  for  the  terminal  engagement  of  the  target  begins 
with  the  missile  seeker  acquiring  the  target.  The  geometry  of  the  missile-target 
engagement  is  shown  in  Figure  9. 


Figure  9.     Missile-target  geometry 

In  describing  the  geometry  of  the  missile-target  encounter,  several 
parameters  must  be  calculated.  The  angle  the  missile's  velocity  vector  makes 
with  the  inertial  reference  frame  is  defined  as  y: 


7m  =  tan 


-i 


ym 


(3-15a) 
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The  magnitude  of  the  missile  velocity  vector  Vm  is  given  by  the  relation: 

Vm=VVi+Vym  (3-15D) 

where  Vym  is  the  vertical  component  of  the  missile  velocity  vector  and  Vxm  is 
the  horizontal  component  of  the  missile  velocity  vector,  measured  in  inertial 
coordinates. 

Similarly,  the  magnitude  of  the  target  velocity  vector  and  and  angle  it 
makes  with  the  inertial  reference  is  given  by 


vt=Vvx2t  +  vy2t 


(3- 16a) 
and 


yt  =  tan  l 


vyt 


(3-16b) 


In  computing  the  variables  used  to  calculate  the  line  of  sight  and  its 
derivatives,  it  is  convenient  to  describe  the  relative  position,  velocity,  and 
acceleration  between  the  target  and  the  missile: 

X  =  Xt-Xm  (3-17a) 

Y  =  Yt-Ym  '                          (3-17b) 

Vx  =  Vxt-Vxm  (3-17c) 

Vy=Vyt-Vym  (3-17d) 
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Ax=  Axt-  Axm  (3-17e) 

Ay=Ayt-Aym.  (3-17f) 

The  line  of  sight  between  missile  and  target  is  given  by  G,  and  is 
defined  as: 

G  =  tan  IxJ  (3_18) 

where  X  and  Y  are  the  relative  cartesian  positions  as  defined  by  equations  (3- 
17a)  and  (3- 17b).  The  first  and  second  time  derivatives  of  O,  denoted  a  and  a, 
are  computed  using  the  relations 


,     (xVy-YVx) 


and 


R2  (3-19) 


.     XAy-YAx     2(xVx  +  YVy)(XYy-YVx) 
R2  R4 


In  an  actual  missile-target  engagement,  it  is  very  unlikely  that  the 
target  line  of  sight  rate  and  acceleration  will  be  measured  accurately  enough  in 
order  to  use  6    and  a  directly  in  the  guidance  algorithm.    An  alternative 
method  is  to  calculate  estimates  of  a   and  a,  denoted  by  6  and  a  ,  based  on 
measured  values  of  C.  This  is  accomplished  through  the  use  of  a  Luenberger 
Observer,  as  shown  by  the  signal  flow  graph  in  Figure  10. 
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Figure  10.  Luenberger  Observer  for  estimating  a  and  G 

This  observer  can  be  described  in  state-space  form  by  the  relation 
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(3-21) 


where  R  is  the  range  from  the  missile  to  the  target,  R  is  the  range  rate,  and  ki, 
k2,  and  k3  are  gains  to  be  determined  by  the  designer  in  order  for  the  estimated 
states  to  closely  follow  the  true  states. 
2.  Missile  Guidance  Laws 

The  objective  of  tactical  missile  guidance  is  to  keep  the  line  of  sight 
angle,  a,  between  missile  and  target  constant.  This  is  desirable  because  if  the 
line  of  sight  remains  constant,  while  the  range  from  missile  to  target  decreases, 
an  intercept  will  occur,  as  shown  in  Figure  11.  Two  missile  guidance  laws  will 
be  presented:  classical  proportional  navigation  and  minimum  time  (Bang-Bang) 
control. 
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/o/ 

/o/ 

^     G  represents  line  of  sight  angle 
between  missile  and  target. 

/^y^ 

<rfT  Missile 

Figure  11.  Constant  Line-of-Sight  Intercept 

a.  Proportional  Navigation   Control 

The  method  most  often  used  in  controlling  homing  missiles  is 
known  as  proportional  navigation,  where  the  magnitude  of  the  transverse 
control  force  is  proportional  to  the  rate  of  change  of  the  line  of  sight  (a). 
Lukenbill  [3]  conducted  research  into  proportional  navigation  control  and  his 
derivation  of  the  proportional  navigation  controller  is  given  below. 

Figure  12  depicts  the  basic  proportional  navigation  scheme. 
Assuming  that  the  seeker  head  of  the  missile  follows  the  target,  the  transverse 
acceleration  perpendicular  to  the  line  of  sight  will  equal  the  acceleration  of  the 
R  vector  in  that  direction.    Mathematically,  the  acceleration  of  R  is 

AR=(R  +  coxcoxR)iR+(2coxR  +  cbxR)ie  (3-22) 


where 


R 


missile/target  line  of  sight  vector 
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R        =         closing  rate  along  R 

R        =         acceleration  along  R 

co        =         angular  rate  of  change  of  R  in  inertial  space 

Ar       =         overall  acceleration  of  R. 

At  this  point,  a  missile  acceleration,  Am,  equal  to  the  target 
acceleration,  At,  will  make  the  line  of  sight  parallel  to  its  original  direction.  As 
long  as  R  remains  along  R  (co=0)  a  missile/target  intercept  is  assured.  So,  the 
transverse  acceleration  command  is 

At-Am=cbxR  +  2(coxR)  (3-23) 
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Figure  12.  Vectorial  Proportional  Navigation  Scheme 

Assuming  the  line  of  sight  rate  is  equal  to  the  angular  rate  of  change  of  R  in 
inertia!  space,  equation  (3-23)  now  becomes 


At-Am=Ra  +  2R6 


(3-24) 


where  a  equal  to  co,  and  a  is  equal  to  ci). 

In  the  classical  proportional  navigation  scheme,  the  missile  course 
is  one  in  which  the  rate  of  change  of  the  missile  heading  is  directly  proportional 
to  the  rate  of  rotation  of  the  line  of  sight  vector  from  the  missile  to  the  target. 
As  a  result,  this  course  change  is  intended  to  counteract  the  rotation  of  the  line 
of  sight,  thus  returning  to  a  constant  bearing  course.    The  movement  of  the 
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missile  and  target  cause  the  line  of  sight  to  rotate,  resulting  in  a  differential 
displacement  between  the  missile  and  the  target  perpendicular  to  the  range  line. 
Figure  13  depicts  this  geometry.  The  proportional  navigation  guidance  law 
attempts  to  generate  an  acceleration  command,  Ac,  perpendicular  to  the  line  of 
sight 

Assume  a  gyro  stabilized  seeker  head,  as  in  the  Sidewinder 
missile.  If  there  is  no  torque  applied  to  the  gyro,  the  seeker  will  not  rotate. 
Assuming  the  seeker  tracks  the  target,  the  gyro  angle  will  follow  the  line  of 
sight.  Applying  the  equation  of  motion  for  a  gyro  stabilized  seeker 

L=Icoft  (3-25) 

where 

L         =         applied  torque 
CO        =         spin  angular  velocity 
I  =         moment  of  inertia  of  the  gyroscope 

Q        =         rate  of  precession  of  the  gyroscope. 
Applying  this  to  the  case  when  the  seeker  head  tracks  the  target,  Q.  is  then 
replaced  by  the  rate  at  which  the  gyro  is  torqued  in  space.  This  is  simply  a, 


Figure  13.  Missile  Acceleration  Orientation 
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Figure   14.  Missile  Acceleration  Relationship 

which  is  the  line  of  sight  rate.  Thus,  equation  (3-25)  becomes 

L  =  Icoc.  (3-26) 

This  torque  is  in  turn  applied  to  the  control  surface  of  the  missile  leading  to  the 
relationship 

Am=kL  =  kIcoa  (3-27) 

where  k  is  a  constant  of  proportionality.  Referring  to  Figure  14,  a  relationship  is 
determined  for  Am  in  terms  of  the  rate  of  change  of  the  missile  flight  path  angle, 
ym.  Given  the  missile  velocity  vector  at  some  point  in  time,  Vm(t),  and  suppose 
the  missile  undergoes  an  acceleration,  Am,  during  an  interval  of  time,  dt.  The 
velocity  vector  is  then  displaced  and  is  represented  by  the  vector  Vm(t+dt).  The 
angle  the  vector  is  traversed  is  simply  dym,  the  differential  missile  flight  path 
angle.  For  small  angles  (which  are  guaranteed  by  making  dt  small)  the 
following  relationship  is  obvious: 

Amdt  =  Vmdym.  (3-28) 
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Dividing  equation  (3-28)  by  the  time  interval,  dt,  the  missile  acceleration  is 
defined  as 

A    =V  ^s.  =  V  7  (3-29) 

^m        ym      ,.  vmlm  W  "'J 

at 

Combining  equations  (3-27)  and  (3-29) 

Vmym=kIcod.  (3-30) 


Dividing  through  by  Vm,  the  proportional  navigation  law  becomes 

Ym  = 


^klco^ 


,  Xn  j 


(3-31) 


or 

Ym=Na.  (3-32) 

Equation  (3-32)  represents  the  classical  proportional  navigation  equation  where 

jm      =         rate  of  change  of  the  missile  heading 

c        =         rate  of  change  of  the  line  of  sight 

N        =         proportional  navigation  ratio. 

The  navigation  ratio  determines  the  sensitivity  of  the  missile 
system.  A  high  navigation  ratio  will  lead  to  rather  high  gains  resulting  in  large 
missile  commands  for  small  changes  in  the  line  of  sight  rate.  On  the  other  hand, 
small  values  for  N  will  lead  to  small  missile  commands  for  a  given  d .  Larger 
navigation  ratios  are  preferred  for  head  on  engagements  and  smaller  ones  are 
preferred  for  tail  chase  cases.  For  this  research  the  navigation  ratio  is  taken  to 
be  four. 
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In  tactical  radar  homing  missiles  using  proportional  navigation 
guidance,  the  seeker  provides  an  effective  measurement  of  the  line-of-sight 
rate,  and  a  Doppler  radar  provides  closing  velocity  information.  [4] 

b.     Bang-Bang  Control  Minimizing  Line-of-Sight  Rate 

In  this  control  law,  the  goal  is  the  same  as  in  proportional 
navigation  control:  to  drive  the  line-of-sight  rate  to  zero,  or  in  other  words,  to 
keep  the  line  of  sight  constant. 

Mathematically,  this  control  law  can  be  stated  as. 


u(t)  =  -Nsign 


(         rrlrrh 


G  + 


V 


0\G\ 


(3-33) 


In  order  to  alleviate  chattering,  a  small  linear  region  will  be  introduced,  similar 
to  the  boost  phase  controller: 


u(t)  = 


-Nsign(a(t)  +  — a(t)|c(t)| 


a(t)  +— d(t)|a(t)| 


>e 


--(a(t)  +  -^-a(t)|a(t)|)      -e<a(t)  +  -^-a(t)|a(t)|<  8.(3-34) 
8  V  2N  J  2N 


In  the  next  chapter,  several  simulations  will  be  presented  in  order 
to  illustrate  the  effectiveness  of  the  Bang-Bang  control  algorithms  as  compared 
to  Proportional  Navigation,  in  both  the  boost  phase  and  in  the  terminal  phase  of 
the  missile  flight 

3.     Missile  Dynamics 

The  signal  flow  graphs  for  the  missile  dynamics  are  shown  in  Figure 
15. 
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Figure  15.  Missile  Dynamics 


In  state-space,  these  signal  flow  graphs  can  be  represented  as 
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IV.VERTICAL  LAUNCH  MISSILE   BOOST  PHASE 

SIMULATION 

A.  OVERVIEW 

In  simulating  the  boost  phase  of  the  vertically  launched  missile,  several 
assumptions  can  be  made  to  reduce  the  complexity  of  the  simulation  to  a  level 
that  won't  detract  from  the  concepts  being  presented: 

1 .  The  missile  is  limited  to  Mach  4. 

2.  The  missile  is  limited  to  30  g's  of  acceleration  in  the  transverse  and 
tangential  directions  combined. 

3.  The  speed  of  sound  in  air  is  constant  at  340  m/s. 

4.  The  density  of  air  is  constant  at  1  Kg/m3. 

5.  The  missile  mass  remains  constant. 

6.  The  lift  and  drag  forces  experienced  by  the  missile  can  be  modeled  by 
equations  (3-6)  through  (3-9). 

B.  MISSILE  PARAMETERS 

At  launch,  the  following  parameters  define  the  state  of  the  missile: 
The  initial  missile  parameters  are: 


Lm=4.2  m 

missile  length 

M=225  Kg 

•    missile  mass 

Ltg=2.1  m 

length  from  center  of  gravity  to  tail 

Lpg=0.15  m 

length  from  center  of  gravity  to  center 

of  pressure 

S=0.99  m2 

missile  reference  lift  surface 
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T=66800  N  missile  thrust 

Icg=2000  kg  m2  missile  moment  of  inertia 

p=l .  kg/m3  density  of  air 

Vmax=1460  m/s  max  missile  velocity 

Additionally,  the  launch  point  will  correspond  to  inertial  reference: 
x(0)=0  m  y(0)=0  m 

Vx(0)=0m/s  Vy(0)=0m/s 

C.  CANNISTER  EGRESS 

It  is  not  possible,  nor  is  it  desirable,  to  have  the  missile  commence  its  tip- 
over  maneuver  until  it  has  egressed  from  the  launch  cannister  and  achieved 
sufficient  distance  from  the  launch  platform  to  minimize  hazards  to  personnel. 
Thus,  the  TVC  actuators  will  have  a  zero  degree  deflection  for  the  first  0.7 
seconds  after  launch.  This  will  ensure  sufficient  missile  altitude  before  tip-over 
begins.  At  time=0.7  seconds  the  missile  position  and  velocity  are: 

x(0.7)=0m  y(0.7)=71.8m 

Vx(0.7)=0  m/s  Vy(0.7)=204.5  m/s 

D.  BOOST  SIMULATION  AND  RESULTS 

The  simulation  was  run  with  four  values  of  the  maximum  thrust  vector 
control  (TVC)  angle:  30°,  45°,  60°,  and  90°.  This  was  done  to  illustrate  what 
effect  a  greater  TVC  angle  will  have  on  the  missile  trajectory  and  velocity. 
Figure  16  shows  the  missile  trajectories  for  the  four  cases. 
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Figure  16.     Missile  Boost  Phase  Trajectories 

Analysis  of  the  trajectories  shows  that  as  the  maximum  TVC  angle  is 
increased,  the  missile  will  achieve  horizontal  flight  at  lower  altitudes  than 
smaller  TVC  angles.  However,  at  angles  greater  than  45°,  there  is  a  significant 
reduction  in  the  downrange  distance.  This  may  be  crucial  if  the  missile  is  to  be 
fired  against  a  target  that  is  close  to  the  ship  at  launch.  At  TVC  angles  less 
than  30°,  the  thrust  is  directed  toward  achieving  missile  altitude  to  a  larger 
extent  than  having  the  missile  achieve  horizontal  flight  in  minimum  time.  This 
may  be  advantageous  if  the  missile  if  fired  against  a  high-altitude  target,  but  not 
against  a  sea-skimming  missile.  Another  consequence  of  having  a  large  TVC 
angle  is  that  much  of  the  thrust  is  spent  rotating  the  missile  to  horizontal  at  the 
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expense  of  missile  kinetic  energy,  which  results  in  a  longer  time  for  the  missile 
to  reach  maximum  velocity.  This  is  shown  graphically  in  Figure  17. 


1400 


Figure  17.     Missile  Boost  Phase  Velocity  Magnitude  Profiles 

At  the  two  highest  TVC  angles  (60°  and  90°),  the  first  several  seconds  are 
spent  turning  the  missile  over  while  not  attaining  much  forward  velocity.  As  a 
result  of  a  higher  turning  rate,  higher  angles  of  attack  are  generated,  as  shown 
in  Figure  18.  This  is  not  desirable  because  at  higher  angles  of  attack,  the  drag 
forces  that  act  to  slow  the  missile  down  increase.  Thus,  it  is  advantageous 
from  an  aerodynamic  point  of  view  to  keep  the  angle  of  attack  as  small  as 
possible  in  order  to  minimize  the  drag  force  acting  on  the  missile. 

Figure  19  shows  the  missile  pitch  angle  as  a  function  of  time.  It  is  of 
interest  to  note  that  there  is  not  a  great  difference  between  total  turning  times 
between  the  30°  TVC  angle  and  the  90°  TVC  angle.  This  is  due  to  the  counter- 
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acting  lift  and  drag  forces  encountered  by  the  missile,  during  high  TVC  angles, 
that  are  caused  by  high  angles  of  attack.  This  indicates  that  increasing  the 
maximum  TVC  angle  beyond  45°  may  achieve  a  somewhat  faster  turning  time, 
but  at  the  expense  of  a  great  loss  of  kinetic  energy  and  therefore  velocity.  This 
is  of  peak  importance,  for  without  sufficient  velocity,  a  missile  has  little  or  no 
chance  of  successfully  intercepting  an  incoming  target,  particularly  a 
maneuvering  target. 
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V.   TERMINAL   PHASE   SIMULATIONS 

In  order  to  compare  the  effectiveness  of  minimum-time  'Bang-Bang'  control 
as  compared  to  proportional  navigation  control,  several  missile-target 
geometries  were  simulated.  Each  geometry  was  run  with  four  control  system 
configurations: 

1.  proportional  navigation  using  analytical  values  of  line-of-sight  rate  and 
acceleration,  which  shall  be  referred  to  as  'analytical  proportional  navigation'; 

2.  proportional  navigation  using  estimates  of  line-of-sight  rate  and 
acceleration,  which  shall  be  referred  to  as  'estimated  proportional  navigation'; 

3.  Bang-Bang  control  using  analytical  values  of  line-of-sight  rate  and 
acceleration,  which  shall  be  referred  to  as  'analytical  Bang-Bang'; 

4.  Bang-Bang  control  using  estimates  of  line-of-sight  rate  and  acceleration, 
which  shall  be  referred  to  as  'estimated  Bang-Bang'; 

A.     CASE  ONE 

In  the  first  set  of  simulations,  the  target  will  be  headed  directly  toward  the 
ship,  flying  a  straight  course  with  no  maneuvers.  The  initial  conditions  are: 
Missile: 


x=0; 

z=1000m; 

Vx=680  m/s  (Mach  2) 

Vz=0; 

Target: 

x=7000  m; 

z=30m 

Vx=-850  m/s  (Mach  2.5); 

Vz=0; 
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Figure  20  shows  the  missile-target  geometries  using  true  line-of-sight  rate 
and  acceleration.  Both  proportional  navigation  and  Bang-Bang  control  was  able 
to  hit  the  target.  Note  that  the  trajectory  using  Bang-Bang  control  flew  a 
straighter  course  to  the  intercept  point  than  the  trajectory  using  proportional 
navigation  control.  This  indicates  that  Bang-Bang  control  is  able  to  guide  the 
missile  onto  an  intercept  trajectory  faster  than  the  proportional-navigation 
controller,  which  is  continuously  changing  the  missile  flight  path  in  order  to 
intercept  the  target. 

Figures  21  and  22  show  control  force  as  a  function  of  time  for  the 
proportional  navigation  controller  and  the  Bang-Bang  controller,  respectively. 
The  proportional  navigation  controller  required  continuous  control  in  order  to 
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Figure  20.     Case  One:  Missile/Target  trajectories  using  true  line 

of  sight  rate 
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correct  the  missile  trajectory  enough  to  enable  the  missile  to  intercept  the 
target.  It  does  not  leave  much  reserve  control  to  counter  a  maneuver  away 
from  the  missile  by  the  target.  Conversely,  the  Bang-Bang  controller  applied 
maximum  control  effort  until  the  line  of  sight  rate  became  zero,  then  shut  off. 
This  allows  the  missile  greater  flexibility  in  case  the  target  maneuvers  in  any 
direction.  If  the  linear  zone  were  not  incorporated  into  the  Bang-Bang 
controller,  however,  the  controller  would  have  chattered  between  maximum 
positive  and  negative  control  force  instead  of  shutting  off.  This  would  have 
decreased  the  missile  velocity  unnecessarily  just  when  the  missile  needs 
maximum  kinetic  energy  for  maneuverability. 
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Figure  21.  Case  One:  Control  force  versus  time  using  proportional 
navigation    control  with  true  line  of  sight  rate 

Figure  23  shows  the  trajectories  using  estimates  of  line  of  sight  rate  and 

acceleration.   It  is  apparent  here  that  the  proportional  navigation  controller  was 

not  able  to  react  fast  enough  in  order  to  hit  the  target.   This  is  due  to  the  delay 

encountered  while  the  observer  filter  was  not  matched  up  with  the  analytical 

values,  as  shown  in  Figure  24.  By  the  time  the  estimated  value  of  line  of  sight 

rate  and  acceleration  matched  the  analytical  values,  the  target  was  at  a  point 

where  proportional  navigation  could  not  generate  enough  of  a  control  signal  fast 
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enough  to  enable  target  intercept.  This  illustrates  the  importance  of  being  able 
to  detect,  virtually  instantaneously,  changes  in  target  velocity  and  direction.The 
Bang-Bang  controller  also  experienced  the  same  delay  before  the  proper 
control  was  applied,  as  shown  in  Figure  25,  but  since  it  immediately  applied 
maximum  control,  the  missile  was  still  able  to  intercept  the  target. 
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Figure  22.     Case  One:  Control  force  versus  time  using  Bang-Bang 
control  with  true  line  of  sight  rate 
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Figure  23.  Case  One:  Missile/Target  trajectories  using  estimated 

values  for  line  of  sight  rate 
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Figure  24.     Case  One:  Control  force  versus  time  using  proportional 
navigation  control  with  estimated  values  for  line  of  sight  rate 
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Figure  25.     Case  One:  Control  force  versus  time  using  Bang-Bang 
control  with  estimated  values  for  line  of  sight  rate 

Figure  26  compares  the  trajectories  of  Bang-Bang  control  using  analytical 

and  estimated  values  of  line  of  sight  rate  and  acceleration.     Because  the 

estimated  values  took  approximately  0.25  seconds  to  match  the  analytical 

values,  there  is  a  period  at  the  start  of  the  simulation  where  the  estimated 

trajectory  has  positive  control  applied  to  it,  while  the  analytic  trajectory  has 

negative  control  applied.    After  the  analytic  and  estimated  values  matched  up, 

they  both  were  able  to  drive  the  line  of  sight  rate  to  zero  in  plenty  of  time  to 

intercept  the  target. 
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Figure  26.  Case  One:  Missile/Target  trajectories  using  Bang-Bang 
control,  comparing  true  and  estimated  values  for  line  of  sight  rate 

B.    CASE  TWO 

In  the  second  set  of  simulations  the  target  will  be  headed  directly  toward 
the  ship,  but  will  initiate  a  12  g  vertical  maneuver  at  problem  start.   The  initial 
conditions  are: 
Missile: 


Target: 


x=0; 

z=1000m 

Vx=680  m/s  (Mach  2) 

Vz=0; 

x=7000  m 

z=30m 

Vx=-850  m/s  (Mach  2.5); 

Vz=0. 
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Figure  27  shows  the  simulation  results  of  the  analytical  proportional 
navigation  and  analytical  Bang-Bang  controllers.  The  Bang-Bang  controlled 
missile  was  able  to  adjust  to  the  12  g  vertical  maneuver  of  the  target  in  order  to 
achieve  intercept.  The  proportional  navigation  controlled  missile  was  not  able 
to  keep  up  with  the  high  g  maneuver  of  the  target,  and  missed  the  target  by 
about  80  meters. 

Figure  28  shows  the  control  force  as  a  function  of  time  for  the  analytical 
proportional  navigation  controlled  missile.  Note  that  the  controller  did  not 
saturate  until  well  into  the  simulation.  This  indicates  that  if  more  control  had 
been  applied  sooner,  the  missile  could  have  intercepted  the  target.  This  is 
illustrated  by  Figure  29,  which  shows  the  control  of  the  analytical  Bang-Bang 
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Figure  27.  Case  Two:  Missile/Target  trajectories  using  true  line  of 

sight  rate 
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Figure  28.  Case  Two:  Control  force  versus  time  using  Proportional 
Navigation  control  with  true  line  of  sight  rate 


Figure  29.  Case  Two:  Control  force  versus  time  using  Bang-Bang 
control  with  true  line  of  sight  rate 

controlled  missile.  Neither  controller  went  to  zero  as  the  Bang-Bang  controllers 
of  Case  One  did.  This  is  because  the  target  is  maneuvering,  constantly 
changing  the  line  of  sight  angle. 

The  simulations  using  estimated  values  for  line  of  sight  rate  and 
acceleration  are  shown  in  Figure  30.  As  was  the  case  with  the  analytical 
expressions,  the  Bang-Bang  controller  using  estimates  was  still  able  to  achieve 
intercept,  while  the  proportional  navigation  controller  missed  the  target.  Note 
that  the  miss  distance  for  the  proportional  navigation  controller  in  Figure  30  is 
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less  than  the  miss  distance  achieved  with  the  analytical  expressions.  This  is 
because  the  target  turned  toward  the  missile,  and  the  delay  caused  by  the 
observer  prevented  the  missile  from  beginning  to  dive  after  the  target  as 
quickly,  thus  leaving  less  of  an  altitude  difference  when  control  is  applied.  If  the 
target  were  to  turn  away  from  the  missile,  the  estimates  of  line  of  sight  rate  and 
acceleration  would  be  worse  than  those  obtained  using  analytical  values. 

Figures  31  and  32  are  the  control  versus  time  graphs  for  the  estimated 
proportional  navigation  controller  and  the  estimated  Bang-Bang  controller.  The 
delay  encountered  by  both  of  the  controllers  is  about  the  same  as  those 
encountered  in  Case  One. 
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Figure  30.  Case  Two:  Missile/Target  trajectories  using  estimated 
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Figure  31.  Case  Two:  Control  force  versus  time  using  proportional 
navigation  control  with  estimated  values  for  line  of  sight  rate 
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Figure  32.  Case  Two:  Control  force  versus  time  using  Bang-Bang 
control  with  estimated  values  for  line  of  sight  rate 
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VI.    CONCLUSIONS   AND   RECOMMENDATIONS 

This  work  has  shown  that  the  minimum-time  (Bang-Bang)  controller  is  an 
effective  algorithm  for  missile  control,  in  both  the  boost  and  terminal  phase  of 
tactical  missile  flight.  It  is  particularly  effective  when  the  target  has  a  speed 
advantage  over  the  missile  or  when  the  target  is  maneuvering. 

As  is  the  case  with  virtually  all  tactical  missiles,  the  accuracy  of  the  target's 
measured  position,  velocity,  and  acceleration  vectors  are  paramount  in 
accurately  predicting  the  parameters  with  which  the  missile  is  controlled  with  a 
Bang-Bang  minimum  time  controller. 

Areas  for  future  study  will  include: 

1.  Simulating  me  entire  missile  flight  from  launch  through  intercept.  The 
focus  of  these  studies  will  be  to  analyze  the  effects  of  missile  velocity,  tip-over 
altitude,  and  target  speed  advantage  on  missile  performance. 

2.  Developing  more  complex  models  in  order  to  better  understand  what 
the  effects  of  a  Bang-Bang  controller  will  have  on  a  particular  missile  system. 
Such  studies  will  use  empirical  aerodynamic  and  physical  data  for  a  particular 
missile,  such  as  SM-II  Block  IV  or  Vertically-Launched  Sea  Sparrow,  and  will 
entail  developing  a  5  or  6  DOF  computer  model. 

3.  Investigating  the  effects  of  noise  on  measurement  of  line  of  sight  and 
its  time  derivatives,  to  the  extent  of  the  effect  noise  has  on  a  Bang-Bang 
controller. 

4.  Investigating  the  effects  of  plant  modelling  errors  on  system  sensitivity 
when  using  Bang-Bang  controllers. 
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APPENDIX   1-  BOOST  PHASE  MATLAB   PROGRAM 


%  VLS  BOOST  PHASE  MISSILE  SIMULATION 

%  T  B  MULL 

%  REVISED  22  APRIL  1992 

clear;clg; 

%  Initial  Missile  Parameters 

Delta=[30  45  60  90]; 


forj=l:length(Delta) 

L_m=3.6; 

M=227; 

L_gt=1.8; 

L_gp=0.15; 

S=0.99; 

T=66800; 

Icg=2000; 

rho=1.00; 

Vmax=340*4; 

dt=.02; 

delta_max=Delta(j); 

%     -- 

%  Initial  Conditions 


%  L_m=  missile  length  in  meters 

%  M=  missile  mass  in  kilograms 

%  L_gt=  length  from  center  of  gravity  to  tail 

%  L_gp=  length  from  center  of  gravity  to  center 

%  of  pressure 

%  S  is  lift  surface  in  square  meters 

%  T=  missile  thrust  in  newtons 

%  leg  is  the  missile  moment  of  inertia  (kg  mA2) 

%  rho  is  the  density  of  air  in  kg/mA3 

%  Vmax  is  the  max  missile  velocity  in  m/s 

%  dt  is  simulation  time  step  size 


tfinal=10.00; 


%  Final  simulation  time. 
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kmax=tfinal/dt; 
time=0.0; 

x=0;xdot=0;xddot=0; 
y=0;ydot=0;yddot=0; 


gamma=90.0; 

gammadot=0; 

theta=90.0; 

thetadot=0; 

X_theta=[thetadot;theta] ; 

X_x=[xdot;x]; 

X_y=[ydot;y]; 

A=[0  0;1  0]; 


%  kmax  is  number  of  time  increments 

%  initial  time 

%  initial  missile  states  (in  x) 

%  initial  missile  states  (in  y) 

%  (these  initial  y  states  correspond  to  missile 

%  position  and  velocity  after  clearing  the 

%  cannister) 

%  gamma  is  inertial  to  velocity  vector 

%  in  degrees  (will  convert  to  radians  in  loop) 

%  theta  is  inertial  to  missile  head 


%  A  and  B  matrices  for  determining  velocity  and 
%  position  of  x,y,  and  theta  base  on  acceleration 
%  commands 


B=[1;0]; 

[phi,del]=c2d(A,B,dt); 
alpha=theta-gamma; 
dtr=pi/180; 


%  converts  A,B  to  discrete  for  simulation 

%  alpha  is  Angle-of-attack 

%  converts  degrees  to  radians 
V_m(l)=sqrt(X_x(l,l)A2+X_y(l,l)A2);  %  V_m  is  missile  velocity  magnitude 
F_drag=0; 
F_lift=0; 
Cd=0; 
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fori=l:kmax-l 
%****algorithm  for  determining  rudder  deflection  angle****************** 

%  test  for  cannister  egress(time<0.6) 
iftime(i)<0.7 

delta(i)=0; 
else 

Theta=X_theta(2,i); 

Thetadot=X_theta(  1  ,i); 

N=T*L_gt/Icg*sin(delta_max*dtr); 

test(i)=Theta+(Thetadot)*abs(Thetadot)/(2*N); 

eps=0.5; 

if  abs(test(i))>eps 

delta(i)=delta_max*sign(test(i)); 
else 

delta(i)=delta_max/eps*test(i); 
end; 
end; 

F_Uft(:,i)=0-5*rho*V_m(i)A2*S*0.1*alpha(i); 

Cd(i)=(2*T/(M*rho*VmaxA2*S)); 

F_drag(:,i)=Cd(i)*0.5*rho*V_m(i)A2*S; 
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%F_drag(:,i)=0; 
xddot(i)=-F_Uft(i)/M*sin(gamma(i)*dtr)+TM*cos(gamma(i)*dtr+delta(i)*dtr)... 

-F_drag(:,i)*cos(gamma(i)*dtr); 
yddot(i)=F_lif t(i  )/M  *cos(gamma(i)  *dtr)+T/M  *  sin(gamma(i )  *dtr+ . . . 

delta(i)*dtr)-F_drag(:,i)*sin(gamma(i)*dtr); 
thetaddot(i)=-F_lift(i)*L_gp/Icg*cos(alpha(i)*dtx)-T*L_gt*sin(delta(i)*dtr)/Icg; 

%  - - 

X_theta(:  ,i+l  )=phi*X_theta(:  ,i)+del*thetaddot(i); 

X_x(:,i+l)=phi*X_x(:,i)+del*xddot(i); 

X_y(:,i+l)=phi*X_y(:,i)+del*yddot(i); 

%  — - 

gamma(i+l  )=atan2(X_y(l  ,i+l  ),X_x(l  ,i+l  ))/dtr; 

%  X_y=[ydot;y]; 

alpha(i+ 1  )=X_theta(2,i+ 1  )-gamma(i+l ); 

delta(i+l)=delta(i); 

V_m(i+ 1  )=sqrt(X_x(  1  ,i+l  )A2+X_y  ( 1  ,i+ 1  )A2); 

time(i+l  )=time(i)+dt; 

end;  %  END  OF  LOOP 

91    -^  ^A*  "A?  :A?  'if  -Af  ^A*  "A*  -A*  -A*  -t*  *A*  -A*  •A*  *A*  *A»  »A>  *A*  *A>  *A*  «A>  •!•  *A*  •!>  *1»  -A*  *1»  -X*  -X*  «1»  »1»  *!•  *A>  *i*  *1*  ~if  *1»  ■»!*  *1*  •!»  »1»  *1»  -X>  -X»  -i*  «i»  *1»  *1*  *!•  *1#  •£>  vL»  -X»  «X*  «1*  «A»  -d*  «X»  *1*  -A> 
i—  *T*  ^  *P  *V*  *P  *p  *T*  "P  *r  *v*  *r  *P  *  ^r  •P'  T  T*  I*  T"  ™  T  ^  *^  *P  T"  ^r  T"  T  ■*■  "P  T"  "I*  "I"  "I"  "P  T"  T"  ^*  T*  ■!■  "I*  T*  T*  "•■  I*  T*  *l*  T*  *P  "t*  I"  T*  1*  "I*  *1*  T*  *»■  T^  "P  T* 

evaKt'X.x'.numlstrG^X.x;']); 

eval([,X_y\num2str(j),,=X_y;,]); 

eval([,V_m,,num2str(j),,=V_m;']); 

eval([,alpha,,num2st^(J),,=alpha;,]); 

eval(['X_theta,,num2str(j),'=X_theta;']); 

delta  max 
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vm=V_m(i+l) 

X=X_x(2,i+D 

Y=X^y(2,i+l) 

Alpha=alpha(i+1) 

Theta=X_theta(2,i+l ) 

end; 

Boost  lplot2a 

save  bocstdata 

clg; 

z=length(X_x(2,:)); 
fork=l:length(Delta) 
eval([,X_x=X_x\num2str(k),Y]); 
eval([,X_y=X_y,,num2str(k);;']); 
axis([0  9000  0  2000]); 
plot(X_x(2,:),X_y(2,:));grid; 
hold  on 
if  k==l 
plot(X_x(2,l:2/dt:z),X_y(2,l:2/dt:z),,o') 
plot(X_x(2,z),X_y(2,z),,o') 
elseif  k==2 
plot(X_x(2,l  :2/dt:z),X_y(2,l  :2/dt:z);o') 
plot(X_x(2,z),X_y(2,z),'o*) 
else 
plot(X_x(2,l:2/dt:z),X_y(2,l:2/dt:z),,o') 


52 


plot(X_x(2,z),X_y(2,z),,o') 

end; 
end; 
hold  off 

title('Missile  Boost  Phase-  Bang-Bang  control'); 
xlabeK'Horizontal  Range  (m)'),ylabel(' Altitude  (m)'); 

%text(0.6,0.8,' 30°  TVC  angleVsc') 

pause;clg; 
axis('normal') 

% 

fork=l:length(Delta) 

eval([,V_m=V_m^num2str(k),,;,]); 

plot(time,V_m); 

hold  on; 
end; 
hold  off 

title('V_m  versus  time');grid; 
xlabel(Time  (sec),),ylabel('Velocity  (m/s)'); 

% 

pause;clg; 

for  k=length(Delta):- 1 : 1 

eval(['alpha=alpha,,num2str(k),,;,]); 

plot(time,alpha),grid; 

hold  on; 
end; 


53 


hold  off 

title('Plot  of  Alpha  versus  time') 

xlabel(Time  (sec)'); 

ylabelC Angle  of  Attack  (degrees)'); 

pause;clg; 

% 


axis([0  5  -20  100]); 
fork=l:length(Delta) 

eval([,X_theta=X_theta',num2str(k),';']); 

plot(time,X_theta(2,:  )),grid; 

hold  on; 
end; 
hold  off 

title('Plot  of  Theta  versus  time') 
xlabel(Time  (sec)'); 
ylabel(Theta  (degrees)'); 
axis('normal') 
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APPENDIX  2-  TERMINAL   PHASE   MATLAB   PROGRAM 

%  terminal. m 

%  TERMINAL  MISSILE-TARGET  ENGAGEMENT 

%  TIM  MULL 

%24  APRIL  1992  1130 

%  This  model  does  not  take  drag  or  lift  into  consideration 

%  -uses  two  different  methods  to  compute  sigma  derivatives: 

%   1 .  calculates  sig_dot,  sig_ddot  using  analytical  expressions 

%  2.  uses  an  observer  to  estimate  sighat,sigdothat,sigddothat 

%  3.  controls  using  sigmadothat,sigmaddothat 

%  4.  utilized  bang  bang  on  sigmadot,sigmaddot  with  a  small  linear  region 

%       in  order  to  alleviate  controller  chatter  once  lock  has  been  achieved. 

clear;clg; 

clc; 

disp(['  ']); 

fork=l:4 

disp([*Control=  \num2str(k)]); 

control=k; 

dt=0.01; 

g/q  sp  *  *  sp  55  !p  ^  5p  ♦  ♦  Initial  C'onditions  ^^^^^^^^^^^^^^^^^^^^^^^^h^^^^^h^^^^^^^^ 
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xm=0; 

vxm=2*340; 

axm=0; 

ym=1000; 

vym=0; 

aym=0; 


%  missile  horizontal  position  in  meters 
%  missile  horizontal  velocity  in  m/s 

%  missile  horizontal  acceleration  in  m/sA2 
%  missile  vertical  position  in  meters 
%  missile  vertical  velocity  in  m/s 

%  missile  vertical  acceleration  in  m/sA2 


xt=5000; 

vxt=-2.5*340; 

axt=0; 


yt=30 

vyt=0 

ayt=0 

time=0; 

% 

%observer  stuff 

pl=-100; 

p2=-100; 

p3=-100; 

K=conv([l  -pl],conv([l  -p2],[l  -p3])); 

kl=K(4); 

k2=K(3); 

k3=K(2); 

X=[0;0;0]; 

% 


%  target  horizontal  position  in  meters 
%  target  horizontal  velocity  in  m/s 

%  target  horizontal  acceleration  in  m/sA2 
%  target  vertical  position  in  meters 
%  target  vertical  velocity  in  m/s 

%  target  vertical  acceleration  in  m/sA2 
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fori=l:1200  %loop  1 

%ifrem(i,100)==0 
%clg; 

%subplot(221),plot(xm,ym) 
%subplot(222),plot(time(  1  :length(U)),U) 
%  end; 
gamma(i)=atan2(vym(i),vxm(i));     %  gamma  is  angle  of  missile  velocity 

%  vector  with  respect  to  inertial  ref 
Vm(i)=sqrt(vxm(i)A2+vym(i)A2);     %  Vm  is  missile  velocity  magnitude 

X(i)=xt(i)-xm(i);  %  X  is  relative  horizontal  position 

Y(i)=yt(i)-ym(i);  %  Y  is  relative  vertical  position 

VX(i)=vxt(i)-vxm(i);  %  VX  is  relative  horizontal  velocity 

VY(i)=vyt(i)-vym(i);  %  VY  is  relative  vertical  velocity 

AX(i)=axt(i)-axm(i);  %  AX  is  relative  horizontal  acceleration 

AY(i)=ayt(i)-aym(i);  %  AY  is  relative  vertical  acceleration 

R(i)=sqrt(X(i)A2+Y(i)A2); 

Rdot(i)=sqrt(VX(i)A2+VY(i)A2); 

%  sigma,sigmadot,sigmaddot  are  the  line 

%  of  sight  between  the  missile  and  the 

%  target,  and  its  1st  2  time  derivatives 
c^**********    NOTF     ****************************************** 

%  This  simulation  is  using  analytical  expressions  for  sigma    * 
%  sigmadot,  and  sigmaddot.  Since  the  target's  velocity  and     * 
%  acceleration  cannot  be  measures  directly,  a  more  realistic    * 
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%  approach  is  to  use  an  observer  to  estimate  sigmadothat  and   * 
%  sigmaddothat.  * 

sigma(i)=atan2(Y(i),X(i)); 

sigmadot(i)=(X(i)*VY(i)-Y(i)*VX(i))/R(i)A2; 

sigmaddot(i)=(X(i)*AY(i)-Y(i)*AX(i))/R(i)A2... 

-2*(X(i)*VX(i)+Y(i)*VY(i))*(X(i)*VY(i)-Y(i)*VX(i))/R(i)A4; 

%    3rd  order  observer  to  obtain  sigma.dot.hat,sigma.ddot.hat 
A=[-k3  1  0 

-k2  -2*Rdot(i)/R(i)  1 

-kl  0  0]; 
B=[k3 

k2 

kl]; 
[phi,del]=c2d(A,B,dt); 


X(:  ,i+l  )=phi*X(:  ,i)+del*sigma(i); 

sigmadothat(i)=X(2,i); 

sigmaddothat(i)=X(3,i); 

%  CONTROL  SECTION 

% 

%  PROP  NAV  CONTROL 
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if  control==l  %  if  loop  1 

U(i)=Vm(i)*4*sigmadot(i); 
if  abs(U(i))>15*9.81  %  if  loop  2 

U(i)=15*9.81*sign(U(i)); 
end;  %  end  if  loop  2 

% 

%  PROP  NAV  CONTROL  based  on  sigmadothat 
elseif  control==2  %  if  loop  1 

U(i)=Vm(i)*4*sigmadothat(i); 
if  abs(U(i))>15*9.81  %  if  loop  2 

U(i)=15*9.81*sign(U(i)); 
end;  %  end  if  loop  2 

% 

elseif  control=3  %  if  loop  1 

%  Bang-Bang  control  based  on  sigmadothat,  sigmaddothat 

N=15*9.81;  %  N  is  maximum  normal  acceration  in  Newtons 

test(i)=sigmadothat(i)+(sigmaddothat(i)*abs(sigmaddothat(i)))/(2*N); 

eps=0.001; 

if  abs(test(i))>eps  %  if  loop  3 

U(i)=N*sign(test(i)); 
else 

U(i)=N/eps*test(i); 
end;  %  end  if  loop  3 

% 

elseif  control==4 
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%  Bang-Bang  control  based  on  sigmadot,  sigmaddot 

N=l 5*9.81;  %  N  is  maximum  normal  acceration  in  Newtons 

test(i)=sigmadot(i)+(sigmaddot(i)*abs(sigmaddot(i)))/(2*N); 

eps=0.001; 

if  abs(test(i))>eps 

U(i)=N*sign(test(i)); 
else 

U(i)=N/eps*test(i); 
end; 
end; 

% 

if  rem(i,50)=0 

dispC  working'); 
end; 

%  Let  Xmsl=[xm(i)  State  space  representation  of  missile  states 

%  vxm(i) 

%  ym(i) 

%  vym(i)]; 

Xmsl(:,i)=[xm(i) 

vxm(i) 

ym(i) 

vym(i)]; 

Am=[0  1  0  0 
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0000 
0001 
0  0  0  0]; 
Bm=[0;-sin(sigma(i));0;cos(sigma(i))];  %  control  force  is  perp  to  sigma!! 

[phim,delm]=c2d(Am,Bm,dt);     %  convert  continuous  state-space  to  discrete 


%  update  missile  states 

Xmsl(:,i+1  )=phim*Xmsl(:  ,i)+delm*U(i); 

xm(i+l  )=Xmsl(  1  ,i+ 1 ); 
vxm(i+ 1  )=Xmsl(2,i+l ); 
axm(i+l)=-sin(sigma(i))*U(i); 
ym(i+l)=Xmsl(3,i+l); 
vym(i+ 1  )=Xmsl(4  ,i+ 1 ); 
aym(i+l)=cos(sigma(i))*U(i); 

%  Let  Xtgt=[xt(i) 
%  vxt(i) 

%  yt(i) 

%  vyt(i)]; 

Xtgt(:,i)=[xt(i) 
vxt(i) 


** 
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yt(i) 

vyt(i)]; 

At=[0  1  0  0 

0000 

0001 

0  0  0  0]; 
Bt=[0;axt(i);0;ayt(i)]; 

[phit,delt]=c2d(At,Bt,dt); 

Xtgt(:,i+l)=phit*Xtgt(:,i)+delt*l; 

xt(i+l)=Xtgt(l,i+l); 

vxt(i+l)=Xtgt(2,i+l); 

axt(i+l)=axt(i); 

yt(i+l)=Xtgt(3,i+D; 

vyt(i+l)=Xtgt(4,i+l); 

ayt(i+l)=ayt(i); 

time(i+ 1  )=time(i)+dt; 


%********  check  for  CPA 

X(i+l)=xt(i+l)-xm(i+l); 

Y(i+l)=yt(i+l)-ym(i+l); 


%  X  is  relative  horizontal  position 
%  Y  is  relative  vertical  position 
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R(i+l)=sqrt(X(i+l)A2+Y(i+l)A2); 

%********  subroutine  to  interpolate  the  2  closest  ranges  to  find  CPA 
ifR(i+l)>R(i) 

xms=[ym(i-l:i+l);vym(i-l:i+l);xm(i-l:i+l);vxm(i-l:i+l)]; 

xts=[yt(i-l:i+l);vyt(i-l:i+l);xt(i-l:i+l);vxt(i-l:i+l)]; 

rl  =  interp(xms(:,l:2),xts(:,l:2)); 

r2  =  interp(xms(:,2:3),xts(:,2:3)); 

Rmin=min([rl  r2]); 
break; 
end; 
end; 
j=k; 

eval([,xm',num2str(j),'  =xm;']); 
evalCrym'.nurr^str^),'  =ym;']); 
evalCrx^nurrLZstrtj),'  =xt;']); 
eval(['yt\num2str(j),'  =yt;']); 
eval([,Rmin',num2str(j),*  =Rmin;']); 
eval([,U,,num2str(j),*  =U;']); 
eval(['save  data\num2str(j)]); 
end; 

clg; 
clear 
load  datal 
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save  data5  xml  yml  xtl  ytl 

clear 

load  data2 

xm2=xm; 

ym2=ym; 

xt2=xt; 

yt2=yt; 

save  data5  xml  yml  xtl  ytl  xm2  ym2  xt2  yt2 

clear 

load  data3 

xm3=xm; 

ym3=ym; 

xt3=xt; 

yt3=yt; 

save  data5  xml  yml  xtl  ytl  xm2  ym2  xt2  yt2  xm3  ym3  xt3  yt3 

clear 

load  data4 

xm4=xm; 

ym4=ym; 

xt4=xt; 

yt4=yt; 
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save  data5  xml  yml  xtl  ytl  xm2  ym2  xt2  yt2  xm3  ym3  xt3  yt3  xm4  ym4  xt4 

yt4 

clear 


%- 


%  This  routine  takes  data  computed  in  terminal.m  and  samples  it  every 
%  8  data  point  in  order  to  make  the  graphs  less  memory  intensive 

clear 
load  data5 

xmla=[];ymla=[];xtla=[];ytla=[];xm2a=[];ym2a=[];xm3a=[];ym3a=[]; 
xm4a=[]  ;ym4a=[] ; 

k=min([length(xml),length(xm2),length(xm3),length(xm4)]); 
for  i=l:k 

ifrem(i,8)=0 

xmla=[xmla  xml(i)]; 
ymla=[ymla  yml(i)]; 
xtla=[xtla  xtl(i)]; 
ytla=[ytla  ytl(i)]; 
xm2a=[xm2a  xm2(i)]; 
ym2a=[ym2a  ym2(i)]; 
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xm3a=[xm3a  xm3(i)]; 

ym3a=[ym3a  ym3(i)]: 

xm4a=[xm4a  xm4(i)]; 

ym4a=[ym4a  ym4(i)]; 
end; 
end; 

axis([0  5000  0  1200]); 

plot(xmla,ymla,,-w,,xtla,ytla,,-w,,xm4a,ym4a,,-w,) 
pause; 
plot(xnL2a,ym2a,,-w,,xtla,ytla,,-w',xm3a,ym3a,,-w') 


clg;clear 

load  datal 

ul=U; 

load  data2 

u2=U; 

load  data3 

u3=U; 

load  data4 

u4=U; 

save  control  ul  u2  u3  u4 

clear 

load  control 

dt=0.01; 
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n=min([length(ul),length(u2),length(u3),length(u4)]); 

tfinal=(n-l)*dt; 

time=0.00:dt:tfinal; 

subplot(21 1  ),plot(time,u  1  (1  :n),,-w'),title('control=l '); 

subplot(212),plot(time,u2(l:n),,-w,),title(,control=2'); 

pause;clg; 

subplot(211),plot(time,u3(l:n),,-w,),title(,control=3'); 

subplot(212),plot(time,u4(l:n),,-w,),title('control=4'); 
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